function diff = mean_fun( n, ind, param)
    % calculate numerical mean of variables
    m = linspace(param.m_l,param.m_u,n)';
    G = (q_fun(m,param)-q_fun(param.m_l,param))./(q_fun(param.m_u,param)-q_fun(param.m_l,param));
    g_p = G(2:end)-G(1:end-1);
    if strcmp(ind, 'wage')
        diff = g_p'*wage_fun(m(1:end-1), param);
    elseif strcmp(ind, 'lnwage')
        diff = g_p'*log(wage_fun(m(1:end-1), param));
    elseif strcmp(ind, 'm')
        diff = g_p'*m(1:end-1);
    elseif strcmp(ind, 'delta')
        diff = g_p'*delta_fun(m(1:end-1), param);
    elseif strcmp(ind, 'x')
        diff = (1-u_fun(param)).*g_p'*m(1:end-1).^(1./(1-param.alpha));
    end

end
    